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Abstract 

An extension of the framework of the Finite Integration Technique (FIT) including dynamic and 
adaptive mesh refinement is presented. After recalling the standard formulation of the FIT, the 
proposed mesh adaptation procedure is described. Besides the linear interpolation approach, a 
novel interpolation technique based on specialized spline functions for approximating the discrete 
electromagnetic field solution during mesh adaptation is introduced. The standard FIT on a fixed 
mesh and the new adaptive approach are applied to a simulation test case with known analytical 
solution. The numerical accuracy of the two methods are shown to be comparable. The dynamic 
mesh approach is, however, much more efficient. This is also demonstrated for the full scale 
modeling of the complete RF gun at the Photo Injector Test Facility DESY Zeuthen (PITZ) on a 
single computer. Results of a detailed design study addressing the effects of individual components 
of the gun onto the beam emittance using a fully self-consistent approach are presented. 
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I. INTRODUCTION 



Memory consumption and CPU time represent the main limitations for large-scale elec- 
tromagnetic field computations. This is especially the case for accelerator physics simula- 
tions involving self-consistent charged particle models based on the so-called Particle- In- Cell 

n 

(PIC) method [l[ . Simulations of this type are an indispensable tool for the design and opti- 
mization of particle accelerators since they offer a full insight into the beam dynamics down 
to the particle level. This is especially important for the simulation of low-energy sections of 
an accelerator, where space charge forces heavily influence the beam but also for, e.g., dark 
current or electron cloud simulations. Apart from the solution of Maxwell's equations on 
a discrete grid space, PIC simulations also include the self-consistent solution of the equa- 
tions of motion. Thus, they provide a full description of the accelerator structure including 
space-charge and wakefield effects {2!. 



In this article we address the important problem of numerical efficiency of such beam 
dynamics simulations using the Finite Integration Technique (FIT). The FIT has been suc- 
cessfully applied for the simulation of a wide range of electromagnetic problems in accelerator 
physics jsl, 0] . We propose an extension of this method including dynamic mesh refinement 
in order to locally adjust the spatial grid resolution according to the dynamics of particles 
and fields. This leads to considerable savings in the overall number of computational degrees 
of freedom and, thus, reduces the computational burden in PIC simulations. 



The article is organized as follows. After describing the governing equations in Sec. [TTl a 
brief review of the FIT on static grids is given in Sec. Illlt which also serves for introducing 
the notation. Sec. IIVI describes the dynamic mesh refinement procedures and Sec. |V] ad- 
dresses specific issues of charged particle simulations on dynamically refined meshes. Sec. IVII 
contains two applications. First, an example with known analytical solution is considered 
for investigating accuracy and efficiency of the proposed method. After, results of a detailed 
design study of the PITZ photo injector using self-consistent PIC simulations are presented. 
The achievements are summarized in Sec. IVIII 
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II. FORMULATION OF THE PHYSICAL PROBLEM 



The electromagnetic part of the physical problem considered is described by Maxwell's 
equations. Their integral form reads 

E-'iS=- I ^B-dA, (2) 
D-dA= [ pdV, (3) 

dV Jv 

B-dA= 0, (4) 

dV 

where A and V denote arbitrary surfaces and volumes. The magnetic field strength is 
indicated by H, B denotes the magnetic flux density, E the electric field strength, D the di- 
electric flux density, and J the electric current density. The electric and magnetic quantities 
are related according to the constitutive equations 

D = eE, (5) 
B = fiH. (6) 

For the linear, isotropic, nondispersive materials considered here the permittivity e and the 
permeability yu are material dependent constants. The conservation of charge follows from 
Ampere's law ([T]) and Gauss' law ([3]) in the form of the continuity equation 



^ I pdV+ I J-dA = 0. (7) 

V JdV 



dt 



It reflects that a temporal change of the total charge, Q = Jy p dV, contained in a volume 
is caused only by a flow of electric current into or out of the volume. 

The link between Maxwell's equations and the motion of charged particles is established 
by the Lorentz force 

F = q(^E + ^x B^ , (8) 
on the one hand and Newton's law 

F = ^P. (9) 

on the other hand. Above, q is the charge of the particle, 7 is the relativistic factor and 
u' = JV, with V the velocity of the particle. The mechanical momentum is denoted by p. 



The equations of motion, thus, read 



'7=JL(E^ly.S), (11) 



dt mo V 7 
where u is the position vector and mo the particle mass at rest. 

III. THE FINITE INTEGRATION FRAMEWORK 

For introducing the notation as well as for completeness this section comprises a short 
review of the Finite Integration Technique. For the spatial discretization of Maxwell's equa- 
tions, it utilizes a staggered, dual orthogonal grid doublet consisting of Np primary 
nodes, which covers the domain of interest Q. The staggered grid arrangement is depicted in 
Fig. [U Throughout this article, we employ the FIT on Cartesian grids in three-dimensional 
space. The following discrete integral state variables are introduced 

E-ds, h:= H-ds, 



b := B-dA, d:= iD-dA, (12) 

J A J A 

] := [ J-dA, q:= [ pdV, 

J A Jv 

where c and A belong to the set of edges and faces of the primary grid (black), and c, A and 
V belong to the set of edges, faces and volumes of its dual (gray). The quantities e and h 
represent grid voltages, whereas b, d and j are field fluxes on the faces of the grid. 

Summing up the voltages e and h along the four edges enclosing one face yields a discrete 
but exact representation of Ampere's and Faraday's law, ([T]) and ([2]), in the form 

4 

dt 



J2±h= (13) 



E±s, = -^s, (14) 

for every face of the computational grid. The orientation of the involved voltages ej, hj 
with respect to the orientation of the loop integral in ([1]) and ([2]) determines the summation 
signs (see Fig. [1]). 



The discrete Gauss' laws are obtained in a similar way by summing up the electric and 
magnetic flux variables d and h deflned on all faces of a cell. Fluxes pointing out of the cell 
are counted positive, incoming fluxes negative. This yields the algebraic set of equations 
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(15) 



(16) 



for every cell of Q and Q corresponding to the Gauss' law for electric and magnetic charges, 
respectively. 

The constitutive equations (|5]) and (|6]) translate to 

2 = M,e, (17) 
b = M^h, (18) 
where and provide an appropriate mapping of grid voltages into fluxes 




FIG. 1. Staggered grid doublet in FIT. The primary grid and its associated quantities are depicted 
in black, the dual grid in gray. 



Collecting the Np voltages and fluxes in the vectors e, h, b, d and j allows for stating 

( IT3|) -( IT6|) for the complete computational domain in the compact form 

d « 



Ce 



dt 
d 



Ch= -d+j 

Sd = q, 
Sb = 0. 



(19) 

(20) 
(21) 
(22) 
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The summation signs are gathered together in the matrices C, C and S, S. These matrices 
are, hence, purely topological. They represent a discrete curl and divergence operator of the 
FIT. It can be shown that these operators maintain the properties of their counterparts in 
continuum [sl in the sense that 

SC = 0, SC = 0, 
CS^ = 0, CS^ = 0. 

This allows, e.g., for the derivation of a discrete continuity equation [cf. Eqn. ([7])] in the 
form 

^q + Sj=0. (23) 

Eqns. (IT9l)-(l22l) are called Maxwell grid equations (MGE). Typically, the leap-frog scheme 
is applied for their integration in time. Alternatively, one may apply the longitudinal- 
transverse split operator scheme introduced in for the time integration, which offers 
better numerical dispersion properties. The latter approach is referred to as LT-FIT standing 
for Longitudinal- Transverse-FIT. 



IV. DYNAMIC MESH REFINEMENT 



We will restrict the discussion to conformal mesh refinements, which maintain the con- 
formity and duality of the Cartesian grid doublet. A refined mesh is obtained by a sequence 
of bisections applied to a given mesh cell. The number of consecutive bisections of a cell is 
referred to as its refinement level L (see Fig. [2]). 

When mesh refinement is applied the discrete quantities assigned to the primary and the 
dual grid need to be recomputed according to the modified mesh. In Fig. |3]the refinement of 
a primary grid cell and the subsequent arrangement of electric grid voltages is illustrated in 
a two-dimensional cut view. The preservation of the duality of the staggered grids requires 
adding new nodes to the dual grid, as well as shifting existing nodes. This is shown in 
Fig. m In order to obtain a value for the new or shifted voltages an interpolation procedure 
has to be carried out. Furthermore, it has to be distinguished between the interpolation of 
voltages oriented in parallel to the refinement [e.g. e^(2') in Fig. Et^b)] and those oriented 
perpendicularly [e.g. e^(2) and e'^{2') in Fig. |3]^c)]. In the following, two interpolation 
procedures, linear and high order spline interpolation, will be discussed. 
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Level 



Level 1 



Level 2 




FIG. 2. Representation of the mesh refinement using consecutive bisections and the refinement 
levels. Non-refined grid cells, such as the cells on the left end have refinement level zero. Each 
bisection increases the refinement level by one, which can be represented by means of a binary tree 
structure. The grid step size at level L reduces as 1/2^. 

A. Linear Interpolation of Grid Voltages 

The application of a linear interpolation for the determination of new or shifted voltages 
is straightforward. The interpolation procedure described below refers to the refinement 
scenario shown in Fig. |3] and |H The lengths of the primary and dual grid edges [cf. ([12])] 
are denoted by |c| and |c|, respectively. Discrete field quantities associated with a specific 
node i or i are denoted as, e.g., c{i). Primed quantities are either new or they have to be 
recomputed during the adaptation procedure. 
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FIG. 3. The refinement of the primary grid and the assignment of electric grid voltages to the 
refined edges is depicted in the x — z plane. In (a) the initial situation is shown. In (b) and (c) 
the middle cell was refined, introducing an additional grid line (dashed). Thus, two nodes 2' and 
6' (circles) are inserted. All voltages that are new or require an update of their value are marked 
with a prime. In (b) the position of the new electric grid voltage e'^{2') is shown. The voltage 
e^(2), oriented along the z-coordinate, has to l^e split into the two voltages 6^(2) and e^(2') as 
shown in (c). 
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FIG. 4. The refinement of the dual grid (gray), imposed by the refined primary grid (thin black 
lines) is depicted. The initial situation is shown in (a). In order to preserve the duality of the two 
grids the dual nodes 2 and 6 have to be shifted. Their new position is indicated by 2' and 6' in (b) 
and (c). In addition, the two new nodes 2" and 6" have to be inserted. In (b) the assignment of 
x-oriented magnetic voltages is shown and in (c) the assignment of z-oriented voltages. 
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a. Refinement The values of new electric grid voltages oriented parallel to the refine- 
ment [e.g. e^(2')] are determined by 

<.(2') = 1M±IM, (24) 

The splitting of voltages oriented perpendicularly to the refinement [e.g. 62(2)] is con- 
strained by the condition 

e:(2)+ e:(2')= e.(2). (25) 

This condition imposes the conservation of the total voltage between the nodes 2 and 3. A 
detailed sketch of this refinement is shown in Fig. [51 In order to assign values to the refined 
voltages e'^(2) and e'^(2') the voltage gradient, i.e., the average electric field, along the axis 
is approximated by a central difference using the neighboring voltages 6^(1) and 6^(3). The 
voltages assigned to the refined edges read 

<,2,H.'(2)|.(e,(2,-^.f|Jl). (26) 

,:(,)H.'(2')|.(M2)H-^.f|^). (27) 

where the sampled electric field variables 6^(2) are introduced as 

the 2;-coordinate of the node i is denoted by z{i), and the A-notation denotes a difference 
of the respective quantity, e.g., A2;(3, 1) = 2;(3) — z[l). Since cells are always divided into 
halves, |c'(2)| equals |c'(2')| and the condition (12^ is fulfilled. 

The interpolation of the magnetic voltages is more cumbersome. Besides the insertion of 
new dual grid nodes, also existing nodes have to be shifted in position in order to preserve 
the duality of the staggered grids. Shifting a node, however, implies a modification of the 
grid voltages along the associated edges. 

The interpolation of magnetic voltages oriented in parallel to the refinement is illustrated 
in Fig. 111(b). The voltages h',j.(Q') and h'^(Q") are given by 

Hl.(6') = K(&) - Az{6,6') . (29) 

A2(7, 5) 

h'S") = hS) + A.(6", 6) ■ (30) 

A2;(7, 5) 
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Similarly to the splitting of electric voltages, the splitting of magnetic voltages, oriented 
perpendicularly to the refinement, is constrained by the conservation of the total voltage. 
However, for a complete description of the refinement algorithm, the bisection of at least 
two neighboring cells has to be considered. This is illustrated in Fig. O For this case, the 
conservation condition for the total voltage between the nodes 1 and 4 can be expressed as 



(31) 



The refined magnetic voltages read 



Ki2') = \c'{2') 



K{2) - z{2) 



z{l) + z{2') 



A/i,(4,2) 



Az(4,2) 

K(2") = A^(2, 2') ■ K{2),+Az{2", 2) ■ h,{3), 

K(3') = \c'i3')\-hM, 

K{3") = Az(3, 3') ■ /i,(3) + Az(3", 3) ■ 

h'M = |c'(4)| 

A/i,(4,2 



M4)H-|£(?)±iM-.(4) 



Az(4,2) 

with the sampled magnetic field hz{i) defined as 



Since the following equations hold true 

z(i) + 2(2') z(2>") + z(A) 



z{2)- 



^(4), 



and 



|c(2)| = |c'(2')|+Az(2,2'), 

|c(3)| = Az(2", 2) + |c'(3')| + Az(3, 3'), 

|c(4)|=Az(3",3) + |c'(4)|, 



the Eqns. (l32])-(l36l) fulfill the condition 



(32) 

(33) 
(34) 
(35) 

(36) 



(37) 



(38) 

(39) 
(40) 
(41) 



12 



(1) 



^41) 



(2) 




3 



\ \ 



2^2'^ 3 
2' 2'' 



3 



4 



(3) 



FIG. 5. Interpolation of electric grid voltages oriented perpendicularly to the refinement. The 
coarse grid voltage e z{2) is split into the two voltages e'^{2) and e^(2') such that the total voltage 
in between the nodes 2 and 3 is preserved [cf. Eqn. ([2^ ]. In order to assign values to e^(2) and 
e^(2') the voltage gradient along the axis is evaluated using 6^(1) and 62(8) [cf. Eqns. ([26]) . (p7|) ]. 



The Eqns. (IMI), (EE]), (EZl), (El), ([301), and (l32])-(l36l) define the rules for performing hnear 
interpolations of the state variables e and h. They enable consistent grid refinement based 
on cell bisection within the FIT framework. 

b. Coarsening In the following, the rules for performing grid coarsening are given. We 
refer again to the Figs. [3] and H] and assume the removal of the nodes 2' and 6' and the 
associated dual nodes. Alongside with the removal of the primary nodes 2' and 6' and the 
connecting edge, the electric grid voltage e^(2') is eliminated from the vector of electric 
voltages. No further modification is required for all electric voltages oriented in parallel to 
the refinement. The constraint f l25l) defines the rule for merging electric voltages oriented 
perpendicularly to the refinement 

e,(2) = e:(2) + e',{2'). (42) 

For the determination of the magnetic coarse grid voltage h^^G), oriented parallel to the 
refinement, a linear interpolation of the neighboring voltages allocated on the refined grid 
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FIG. 6. Interpolation of magnetic grid voltages oriented perpendicularly to the refinement. The 
voltages hz{2), hz{3) and hz{i), assigned to coarse grid edges, are split and distributed such that 
the total magnetic voltage in between the nodes 2 and 4 is preserved [cf. Eqn. ()3ip ]. The voltages 
on the refined grid are determined according to the [cf. Eqns. (|32p -(j36 p ]. 

is performed 

E.(6) = (43) 

The magnetic coarse grid voltages, oriented perpendicularly to the refinement, are ob- 
tained by the summation of fractional voltages given on the fine grid (see Fig. [6]). The 
merging of the fine grid magnetic voltages has to obey the constraint fl3Tl) . The coarse grid 
voltages are given by 

hz{2) = K(2') + Az{l2') . K(2''), (44) 
hz{3) = Az{2",2)-K{2") 

+ Ki3') + Az{3,3')-K(3"), 
= Az(3~',3)/i;(3") + h'M. (46) 



(45) 



Since the following holds true 

\c{2")\ = Az(2, 2') + Az{2",2), (47) 
|c(3")| = Az(3, 3') + Az(3", 3), (48) 

the condition (|3T1) is fulfilled also for the case of grid coarsening. 
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Eqns. (??)- (H6|) describe the required modifications of tlie FIT field quantities during 
the coarsening process of the computational grid. In combination with the refinement rules 
given above, they provide a consistent framework for performing time-adaptive conformal 
grid refinement within the FIT. Refinement schemes other than by cell bisecting are generally 
possible. However, these are associated with a larger complexity in mesh administration and 
will not be considered in this work. 



B. Spline Interpolation 

A linear interpolation of the discrete field quantities is easy to implement and fast in code 
execution. The interpolated quantities are second order accurate, which is in agreement with 
the theoretical order of accuracy of the FIT. However, applying higher order interpolating 
functions may be beneficial for obtaining a smoother representation of high-frequency fields 
on the refined mesh. 

Polynomials offer a possibility for performing higher order interpolations. However, poly- 
nomials of high degrees tend to exhibit an oscillatory behavior and possibly large overshoots, 
known as Runge's phenomenon Q|. Commonly, spline functions are employed in order to 
avoid this behavior. 

A spline is a piecewise defined polynomial function of order P which is P — 1 times 
continuously differentiable, i.e., S G C^~^ 9|. The construction rules determine the specific 
type of spline. If the spline S is demanded to pass exactly through the given data points, 
and to be twice continuously differentiable (i.e. S G C^) with the second derivative equal to 
zero on every interval boundary, the so-called natural cubic spline (C-spline) is obtained. 
For its determination, a tridiagonal system of equations has to be solved. During a time- 
domain simulation adopting time- adaptive grid refinement, hundreds to thousands of grid 
adaptations are performed, each involving thousands of cells. Thus, the solution of a system 
of equations is computationally very expensive. In addition, the natural cubic spline may 
still exhibit overshoots (see Fig. [7]). 

In order to further mitigate Runge's phenomenon, the conditions on the continuous dif- 
ferentiability of the spline have to be reduced. A spline S of order P which is at most P — 2 



times continuously differentiable is called a broken spline or subspline 



is a cubic subspline which is an element of the functional space [10 1 
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9|. The Akima spline 
Cubic (sub-) splines 



can conveniently be characterized by the triple of values [x{i), f{i), s{i)] for each data point 
i, where x{i) is the coordinate of the data point, f{i) its value, and s{i) its slope. Imposing 
the conditions 



S[x{-l)] = f{t-l), %(z)] = /(z), 

S'[x{i - 1)] = s{i - 1), S'[x{i)] = s{i), 

at every interval boundary, uniquely describes the broken cubic spline function 

S{x) = ao + ai[x — x{i — 1)] + a2[x — x{i — 1)]^ + a^lx — x{i — 1)]^, (49) 

with the coefficients 

ao = f{i - 1), 
ai = s{i - 1), 

«2 = (^^/(^' ^ - 1) - ^[^(0 + - 1)], 

a3 = ^A/(M-l)+(^A.(M-l). 

However, the slopes s{i) are, in general, unknown. The Akima procedure makes use of a 
local heuristic method for the estimation of the slopes. It involves the data point i and two 
neighboring points on each side. First, the piecewise gradients g, given by 

'^'^'■= xU + l)-xU) ' with j=z-2,z-l,z,z + l, 
are evaluated. These are weighted with the factors w 

wii-1) := \gii + 1) - g{z)\, 
wiz):= \gii-l)-git-2)l 

yielding the estimated slope at data point i 

^^^^ w{^-l)9{^- l) + w{^)g{^)^ (50) 
w{i-l) + w(i) ^ ' 

The fundamental idea of the heuristic can be summarized as: the larger the difference 
between the two gradients on one side of point i, the larger is the weighting factor applied 
to the gradient on the other side of this point. This is intended to minimize overshooting 
and oscillatory behavior. Note that, Akima's heuristic is closely related to the idea of slope 
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FIG. 7. Comparison of spline interpolations. The given data points are interpolated using linear 
interpolation, a natural cubic spline (C-spline) and an Akima spline. The C-spline exhibits an 
oscillatory behavior and strong overshooting. Up to the data point 10, Akima's original slope 
estimation was employed [cf. Eqn. (I50p ]. From this point on the minmod slope limiter [cf. Eqn. (I5ip ] 
was applied, which effectively avoids overshooting. 

limiters. These are commonly used in Finite Volume methods for reconstructing a piecewise 
continuous solution on the grid (cf. [ll|). Using, e.g., the minmod limiter the estimated 
slope reads 



s(i 



minmod[g'(z), g{i + 1)] 



g{i), {ii I \g{i)\ < + + >0}, 

g{i + l), {ii \\g{i)\>\g{i + 1)1 gii)gii + l)>0}, (51) 
0, {Vi\g{i}g{i + 1)<0}. 

Out of its two arguments, the minmod operator chooses the smaller one if they are equally 
signed and zero otherwise. An overview of slope limiting techniques is found in 12 1. 

In Fig. [7] the interpolation of a set of data points using linear interpolation, the C-spline 
and the Akima spline is illustrated. The C-spline shows an oscillatory behavior and strong 
overshooting. For the determination of the Akima spline up to data point 10, the slope was 
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1 2 3 3' 4 4' 5 5' 6 7 



FIG. 8. Setup and evaluation of the interpolating spline. An exemplary arrangement of grid nodes 
along the z-coordinate is shown in (a) and (b). In (b) the refined area of the grid has moved by a 
distance of one grid cell into the z-direction. Hence, the grid has to be coarsened around node 2 
and refined around node 5. The interpolating spline is set up, using the position of the nodes and 
field data in the refined area and some neighboring cells (black dots). An evaluation of the spline 
function at the position of new or shifted nodes (circles), yields the interpolated field values. 

calculated using Akima's slope definition f l50|) . For the points 11 to 16 the minmod limiter 
was applied, which effectively avoids any overshooting. 

The actual spline interpolation procedure of the discrete field quantities is illustrated in 
Fig. [HI In (a) and (b) the primary and dual grid nodes along the z-coordinate are indicated 
in black and gray. In (b) the refined area has moved by a distance of one grid cell into the 
z-direction. Therefore, the grid topology around the cells 2 (coarsening) and 5 (refinement) 
has to be modified. The spline is set up using the position of the grid nodes and the allocated 
field data in the refined area and some neighboring cells. An evaluation of the spline function 
at the position of new or shifted nodes yields the interpolated field values. 

The Akima and the slope limited splines are viable choices for performing interpolations 
within the computational grid. They are formally third order accurate and can be fully 
constructed from local information. Thus, they do not require the solution of a system of 

18 



equations, which makes their implementation very efficient for the purpose of adaptive grid 
refinement. 



V. CHARGED PARTICLE SIMULATIONS ON ADAPTIVE GRIDS 



Simulations of the dynamics of charged particles using the FIT are performed routinely 



and have ffist been reported in 1988 in 



15|. For self-consistent simulations the particle- 



in-cell (PIC) algorithm is applied [16|. It employs macro particles, which may adopt any 
position from the space-continuous domain of interest. Macro particles carry the charge 
and mass of about lO^ '^ individual particles. Due to flTT]) the trajectory of a particle in 
the presence of electromagnetic fields depends only on the ratio of its charge and rest mass, 
making this a valid simplification. 

The PIC algorithm consists of three steps. First, the electromagnetic field at the space- 
continuous position of each macro particle has to be obtained from the space-discrete FIT 
quantities. This is done by means of a trilinear interpolation. Next, the equations of motion 
ffTOj) . f|TT]) are integrated using the algorithm described in [l3|. Finally, the convective cur- 



rents are calculated from the position increment of all particles. For this last step we adapted 



the cloud-in-cell (CIC) approach [18|, ll9[ to work with time-adaptive mesh refinement. The 
particles are modeled as a uniformly charged volume of finite extent, i.e., a charged could. 

The extent of the cloud is usually chosen to coincide with the size of the cells. However, if 
a nonequidistant grid is employed or in the case of adaptive grid refinement, the nonuniform 
sizes of the cells demand for a modification of the algorithm. There are two modification 
options: ffist, a constant size of the particle cloud is chosen independently from the grid 
cell size or, second, the size of the cloud is adapted in order to match the sizes of the 
involved cells. Depending on the local degree of grid refinement, the ffist option can largely 
increase the number of cells affected by one cloud. Besides the coding efforts coming along 
with this non-local operation, the deposit of fractional currents to a large number of grid 
points increases the computational load. Therefore, the second option has been pursued. 
The adaptation of the cloud in dependence of the grid cell size is illustrated in Fig. [9] for a 
one-dimensional example. 

The modified CIC approach maintains its charge conserving property, meaning that it 
fulfills the discrete continuity equation (|23|1 in a cell-wise manner. 
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FIG. 9. Illustration of the adaptation of the particle cloud size in one dimension. The size of the 
particle cloud relates to the extent of the covered cells on the dual grid. If the particle is centered 
within a dual grid cell (cases 1, 3, 5, 7 and 9) the sizes of the cloud and the dual volume coincide. 
Otherwise, the cloud has to be adapted asymmetrically around the particle position. For the cases 
2 and 4 the charge density within the part of the cloud within the dual volumes 2 or 3' has to be 
increased as a consequence of the diminished size. Contrarily, for the cases 6 and 8, the charge 
density of the cloud in the dual volumes 3" and 4 has to be decreased. 



VI. APPLICATIONS 



In this section, we present results of the application of the dynamic mesh adaptation 
algorithm introduced above to the simulation of two setups. 
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A. Bunch drift in a pipe 



We consider a hollow, perfectly conducting pipe, which is closed by a perfectly conducting 
plate at one end. The other end is left open. A bunch of charged particles is emitted from 
the center of the end plate and travels along the pipe at a constant velocity of 0.9 Cq, where 
Co is the speed of light in vacuum. The settings of the computational model are specified in 



Table [H In 20| the analytical solution of this problem for a semi-infinite pipe is given. 

TABLE I. Settings of the benchmark example 



Pipe length 


Pipe radius 


RMS bunch radius 


RMS bunch length 


Particle velocity 


120 mm 


40 mm 


5 mm 


3 mm 


0.9 Co 



In the simulations the axis of the cylinder is aligned with the ^-coordinate. In the 
Figs. [TOT a) and lTOT c). the analytical solution of the longitudinal electric field along the axis 
is depicted by a black curve. The solutions obtained with the FIT and the LT-FIT are 
shown in red and green, respectively. For the results shown in Fig. ITOlfa). an equidistant 
grid and the respective maximum time step was applied. In Fig. [TOT b) the errors, given by 
the pointwise Euclidean distance of the analytical and the respective numerical solution, are 
plotted. The results obtained using dynamic mesh refinement (L = 3) in the bunch region 
are shown in Fig. [TOT c) and the errors in Fig. [TOT d). 

A series of simulations of this test setup using various mesh and mesh refinement settings 
has been carried out. The results are presented in the Tables HTl (FIT) and IIIII (LT-FIT). 
Besides the mesh settings, they list the number of degrees of freedom (DoF) applied as well 
as the computation time in seconds. The relative L^-error and the total variation (TV) of 
the longitudinal electric field along the cylinder axis are given as measures for the quality 
of the numerical solutions. The relative L^-error is computed as 

llE.lh ' 

with = [^Ez{z{l)), ..Ezi^z^Nz)))"^ and the longitudinal component of the analytical 
solution of the electric field evaluated at the grid point positions z{izY 



The total variation is a measure for the smoothness of a function 



211 ]. It is defined as 



TViE{z)) := / \E'{z)\dz. (52) 
^(1) 
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TABLE II. Results of FIT simulations: the upper half shows results for static, equidistant grids, 
in the lower half the results obtained on time-adaptive grids using Akima splines are displayed. 



L 






Ax / mm 


Azmin / mm 


DoF / le6 


Time / sec 




TV 


static grid 





135 


210 


0.59 


0.57 


22.96 


2066 


0.041 


1.00 


time-adaptively refined grid 


1 


135 


105 


0.59 


0.57 


13.00 


1225 


0.044 


1.42 


2 


135 


52 


0.59 


0.58 


9.00 


940 


0.040 


1.93 


3 


135 


27 


0.59 


0.56 


8.50 


925 


0.041 


1.99 


4 


135 


13 


0.59 


0.58 


11.50 


1025 


0.039 


1.55 



TABLE III. Results of LT-FIT simulations: the upper half shows results for static, equidistant grids, 
in the lower half the results obtained on time-adaptive grids using Akima splines are displayed. 



L 






Ax / mm 


Azmin / mm 


DoF / le6 


Time / sec 




TV 


static grid 





135 


210 


0.59 


0.57 


22.96 


3543 


0.038 


0.63 


time-adaptively refined grid 


1 


135 


105 


0.59 


0.57 


13.00 


2008 


0.045 


1.39 


2 


135 


52 


0.59 


0.58 


9.00 


1443 


0.041 


1.65 


3 


135 


27 


0.59 


0.56 


8.50 


1447 


0.040 


1.80 


4 


135 


13 


0.59 


0.58 


11.50 


1791 


0.037 


1.40 



In 22[ it is shown that the TV of a function equals the sum of local minima and maxima, 
where the values at the integration endpoints count once and all other extrema count twice. 
The TV is, therefore, a direct measure for oscillations and their amplitudes. For discrete 
solutions it is computed as 



TV{e,) = J2 |e.(^. + l) -e, 

iz=l 



(53) 



The standard FIT is a very well established method for performing PIC simulations. 
All TV values are, hence, normalized to the value obtained with the FIT on the finest. 
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FIG. 10. Longitudinal component of the electric field, excited by a Gaussian bunch traveling in 
a pipe. The black curve in (a) and (c) indicates the analytical solution for the setup described in 
Table HI All other curves shows results obtained by the numerical simulation of the problem. In 
the top row, a static, equidistant computational grid was applied. The results are shown in (a), 
their errors in (b). In the bottom row, time-adaptive grid refinement (L = 3) was applied. In (c) 
the results for this case are shown and the errors in (d). 

nonadaptive grid, which has been employed. For comparabihty, the settings were chosen 
such that all simulations finish within less than one hour. 
A comparison of the results shows that: 

- The relative errors of the FIT and the LT-FIT are similar, however, the TV of the 
LT-FIT solution is significantly lower, indicating reduced oscillations due to better 
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numerical dispersion properties. 

- For identical numbers of DoF the LT-FIT simulations take longer because of the higher 
computational costs for the time integration. 

- For identical mesh resolutions in the bunch area the adaptive simulations yield errors 
similar to those obtained on nonadaptive meshes. This follows from a comparison 
of the result of the adaptive simulations with the respective nonadaptive simulation. 
Hence, the accuracy of the solutions for this example is mainly determined by the 
resolution in the bunch region. 



B. Self-consistent simulation of the PITZ RF gun 

In order to drive a free-electron-laser (FEL) operating by the self-amplified spontaneous- 
emission (SASE) principle, highly charged electron beams of high brightness are re quir ed 23|. 
The aim of the PITZ project (Photo Injector Test Facility at DESY Zeuthen) [2J] is the 
development and testing of an injector capable of delivering such high quality beams. The 
injectors for the Free Electron Las er in Hamburg (FLASH) 25| and the future European 
X-Ray Laser Project XFEL 26|] are under development at PITZ. 

The layout of the radio-frequency (RF) gun of the injector is shown in Fig. [TTl The 
emitted bunch is accelerated in a 1.5-cell L-band cavity providing for an accelerating gradient 
of 42 MV/m at an operational frequency of 1.3 GHz. A focusing technique proposed in [^^j is 
applied in order to compensate for the correlated part of the space charge induced emittance 
growth. The initial gun layout was designed such that a minimum of the transverse emittance 
is expected at a distance of approximately 1.6 m downstream of the cathode. At this position 
a RF cavity is installed in order to accelerate the electrons to relativistic energies. The main 
design parameters are listed in Table IIVI 



Many numerical studies of the injector using P 



last years. Results have been published, e.g., in [28 



C codes have been performed over the 
-|32|. However, these simulations either 



assumed a rotationally symmetric geometry and made use of a so-called 2.5-dimensional 



approach as implemented in the MAFIA TS2 module 33|, or otherwise very short distances 



in the cm-range of the full three-dimensional model were simulated. Results of parallelized 
simulations of the full three-dimensional model up to one meter downstream of the cathode 
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TABLE IV. PITZ design parameters 



Parameter Design value 

Bunch charge 1 nC 

Transverse laser profile Hat-profile 

Laser spot radius 1 mm 

Longitudinal laser profile Flat-top 

Laser pulse duration 22 ps 

Rise/Fall time 2 ps 

Accelerating gradient 42 MV/m 

Transverse emittance 1 mm mrad 



have been presented in j^, |3^. However, the exact position of the transverse emittance 
minimum is required for the optimal placement of the accelerating cavity. Also, the value 
of the emittance at this position is of strong interest in order to check whether the design 
value is met. 

In the context of this work, the gun was simulated up to a distance of two meters down- 
stream of the cathode. First, settings of the simulation parameter such as the grid resolution 
and the number of macro particles have been determined. Then, a design study was per- 
formed, which addresses the effects of individual injector elements on the beam quality. 

Previous investigations have shown that the length of the simulated bunches critically 



depends on the longitudinal grid step size 



35|. In 



it was stated that a grid resolution 



of 20 /xm longitudinally is required in order to obtain accurate results. In the Fig. [TW a) the 
results of a parameter study addressing the relation of the longitudinal grid step size and the 
computed bunch length are presented. The initial 2.5 cm of the gun were simulated using 
an equidistant grid of 2.5 /im to 80 fim step size. The computed RMS bunch length at 
z = 2.5 cm varies from 2.24 mm to 2.35 mm. Changing the step size from 5 yum to 2.5 fim 
results in a variation of by approximately 1 fim. The computed bunch lengths for a step 
size of 2.5 /im and 10 /im differ by approximately 10 /im. A longitudinal grid step of 10 /xm 
is considered to be a reasonable balance of computational costs and accuracy. 

After emission the electrons have a very low energy of ~ 5 eV and space charge forces have 
a strong influence. Since the particles gain energy quickly, it is sufficient to apply the smallest 
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grid step size only in the immediate vicinity of the cathode. Hence, an additional static grid 
refinement is applied in this area. In the Fig. [T2r b) results for this approach are shown. 
A grid with a uniform step size of 80 fim is statically refined within the first centimeter 
from the cathode. Refinement levels from one to four have been applied, resulting in step 
sizes of 40, 20, 10, and 5 fim within the refined region. The results are, except for minor 
discrepancies, identical to those obtained with an equidistant grid of the same minimum 
step size. 

For the simulation of the full structure, a combination of static and dynamic mesh refine- 
ment was applied providing for a longitudinal resolution of 10 /im for the first cm from the 
cathode and 80 yum thereafter. The requirements on the transverse grid step size are less 
demanding. The width of the simulated bunches using a step size of 2.5 fim and 80 fim differ 
by less than 5 /im. Deviations in the |Um-range correspond to some per mille of the bunch 
width. Choosing a minimum transverse step size of 80 fim and a refinement level of five, 
we obtained an adaptive computational mesh consisting of an average of approximately 82 
million cells. A nonadaptive mesh providing for the same minimum resolution would consist 
of approximately 1.7 billion cells. The bunches in the simulations comprised approximately 
0.5 million macro particles. 

The results given in the following were obtained with the LT-FIT method. In Fig. [131 the 
evolution of the horizontal root-mean-square (RMS) bunch width (a), the RMS bunch 
length az (b), the average particle energy (c), and the horizontal emittance (d) for the 
model shown in Fig. [11] are plotted. The projected emittances and Sy are given by 
Eu = {crl ■ CTp^ — G {x,y}, where 0"^^ is the momentum variance and cru,pu is 

the covariance of position and momentum. 

We performed a design study of the injector in order to identify the individual effects of 
the diagnostics section, the laser mirror, and the shutter valve on the beam quality in terms 
of emittance growth. The respective models are shown in the Figs. [TW a)-(d). For each 
model, an additional element is added. The colors of the curves in Fig. [TW e)-(f) correspond 
to the outline color of the respective model. 

For the model given in Fig. [T^ a). the horizontal and vertical emittances, and Ey, 
are identical. This is expected since the geometry is rotationally symmetric for all parts 
close to the beam. The RF input coupler is obviously not rotationally symmetric but this 
symmetry violation is hidden by the coupling antenna. In Fig. [MT b). the large opening on 
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the bottom side of the doublecross introduces an asymmetry of the geometry in horizontal 
and vertical direction. While Sx is almost unchanged, the value of Sy increases by approxi- 
mately 0.04 mm mrad. The laser mirror, added in the horizontal plane, introduces another 
asymmetry in Fig. ITW c). An influence on the transverse emittances can be identified in the 
results but the actual emittance growth at the position of the minimum is small. Finally, 
the shutter valve depicted in Fig. [T^ d) results in significant changes of the transverse emit- 
tances. While the horizontal emittance actually decreases by approximately 0.1 mm mrad, 
the vertical emittance increases by about 0.18 mm mrad in comparison to the structure of 
Fig. [T^ a). Hence, the shutter valve causes a distinct asymmetry in the transverse beam 
dynamics. In the meantime, a shutter valve with RF shieldings has been installed. 



resonant cavity p.^ j^put coupler 




FIG. 11. A CAD model of the injector section of the Free-electron LASer in Hamburg (FLASH) 
is shown in cut view. The depicted part corresponds to approximately 60 cm of the injector. In 
total 2 m have been modeled. However, the non-depicted part consists only of the beam pipe. A 
short laser pulse is directed onto the photo cathode, where electrons are emitted. The bunch of 
electrons is accelerated by a high-frequency electromagnetic field, which is excited in the resonant 
cavity. It propagates along the electron path through the coupling antenna and passes the shutter 
valve and the laser mirror, which is inserted through one arm of the diagnostics section. The other 
arms are used for inserting measurement devices. The opening on its bottom side is connected to 
a vacuum pump. 
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(a) (b) 




0.5 1.0 1.5 2.0 2.5 0.5 1.0 1.5 2.0 2.5 



z I cm z I cm 

FIG. 12. Computed RMS bunch length vs. longitudinal grid step size. In (a) equidistant grids 
with different grid step sizes were used. In (b) the equidistant grid with a step size of 80 iim. was 
statically refined within the first cm from the cathode. The step sizes given in the legend are in 
units of [xm. 



VII. CONCLUSIONS 



In this paper, we extended the framework of Finite Integration Technique to include 
dynamic mesh refinement. This offers computability on workstations for a class of multi- 
scale problems, which before was accessible to massively parallelized computations only. In 
particular, our approach enabled us to perform fully self-consistent simulations of the PITZ 
RF gun up to two meters downstream from the cathode. We elaborated on the details 
concerning the mesh refinement and coarsening procedures and presented a novel type of 
sub-spline interpolation, which is entirely devoid of overshooting. Finally, we performed a 
design study, which identifies the emittance growth due to individual parts of the gun. 
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FIG. 13. The plots show the evolution of the RMS bunch width (a), RMS bunch length (b), 
average particle energy (c), and horizontal emittance (d) along the longitudinal coordinate of the 
PITZ RF gun. 

analytical solution to the example used in Sec. IVI A[ 
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FIG. 14. Models and results of the PITZ RF gun design study. In (a) to (d) the CAD models 
of the simulated structures are shown. The horizontal and vertical emittances around the position 
of their minimum are plotted in (e) and (f) respectively. The colors of the curves correspond to 
the outline color of the models above. While £x and £y are similar for the case (a), they show a 
distinct asymmetry for the model depicted in (d). 
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